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SUMMARY 


A vectorized, semi-implicit code is developed for the 
solution of the time-dependent, three-dimensional equations of 
motion in plane Poiseuille flow by the large-eddy simulation 
technique. The code is tested by comparing results with those 
obtained from the solutions of the Orr-Sommerfeld equation. 
Comparisons indicate that finite-differences employed along 
the cross-stream direction act as an implicit filter. This 
removes the necessity of explicit filtering along this 
direction (where a nonhomogeneous mesh is used) for the simu- 
lation of laminar flow transition which includes nonlinear 
effects . 


1. INTRODUCTION 


This report is a summary of work on the computation of 
laminar flow control {LFC) by suction and laminar flow transi- 
tion. We consider the time— dependent , three-dimensional equa- 
tions of motion for incompressible, plane Poiseuille flow and 
solve these equations numerically by the large-eddy simulation 
technique. The main features of the solution procedure 
involve the use of the pseudo-spectral technique to evaluate 
spatial derivatives along homogeneous directions and finite 
differences along the nonhomogeneous direction, solving the 
continuity equation directly and transformation of the finite 
difference equations into wave-number space to reduce the 
partial differential equation set to a set of ordinary differ- 
ential equations which are solved by block tridiagonal matrix 
inversion. It should be noted that although the present 
approach makes use of the large-eddy simulation technique which 
enables it to be valid for both laminar and turbulent flows, 
during this contract period we have performed calculations only 
for laminar flow stability for code testing. The incorporation 
of a subgrid scale turbulence model, which can be conveniently 
done by simply changing one subroutine, is necessary for turbu- 
lent flow computations. 

It is generally accepted that transition of laminar flow 
to turbulent flow involves a number of intermediate stages. The 
initial stage is characterized by the existence of two-dimen- 
sional Tollmien-Schlichting v/aves which become weakly three- 
dimensional in the second stage. This is followed by the third 
stage in which the Tollmien-Schlichting waves become fully 
three-dimensional. In the fourth stage high shear-layers with 
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inflexional velocity profiles are developed and these shear- 
layers produce secondary instabilities in the fifth stage. 
Finally, the last stage is the breakdown of three-dimensional 
wave pattern into turbulence. 

Mechanisms and phenomena that follow the amplification of 
the Tollmien-Schlichting waves and that finally lead to the 
fully developed turbulent flow cannot be predicted by linear 
theory. The theoretical understanding and prediction methods 
for this increasingly three-dimensional flow in which nonlinear 
effects play a dominant role do not show close agreement with 
experimental observations. Hence, although the behavior of 
small amplitude disturbances on laminar flow stability is well 
understood, there is as yet no close agreement between theoret- 
ical and experimental studies of the behavior of finite-ampli- 
tude disturbances and their effect on transition. 

This work involves a complete simulation of laminar flow 
transition to turbulence, including nonlinear effects and as 
reported here involves the following tasks: 

(1) Vectorization of the channel flow solution procedure of 
reference 3 for the CYBER-203 computer at NASA/Langley Research 
Center . 

(2) Development and formulation of various initial conditions 
required for the calculation of the transition and LFC. 

(3) Preliminary calculations for laminar flow stability 
for code testing purposes. 

(4) Definition and implementation of required wall-boundary 
conditions for LFC by suction. 

(5) An initial study of subgrid scale model and length 
scale . 
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All the calculations reported herein were performed on 
a 32x33x32 mesh. The objective of these calculations have been 
to test the code, e.g. with regard to speed and efficiency, as 
well as accuracy and stability for various boundary conditions. 
Calculations were performed with arbitrary periodic initial 
conditions that satisfy no-slip conditions and with initial 
conditions that are the eigensolutions of the Orr-Sommerfeld 
equation. These are reported in section 3, whereas section 2 
contains a brief description of the solution procedure. 

2. SOLUTION PROCEDURE 

In large-eddy simulation of turbulent flows, each flow 
variable is decomposed into a large-scale component and a 
residual-field component. This can be obtained by filtering 
the flow variables, u^ and p (where are the velocity 
components along x^ and p is the pressure) such that 

u^(x,t) = u^ + u'(x,t) (1) 


and 


u^ (x,t) 


G (x-x ' ) u . (x ' , t) dx ' 
J ~ ~ 1 

D 


( 2 ) 


In these expressions u^ is the residual-field component, 
u^ is the large-scale (filtered) component, and G(x,x') is the 
Gaussian filter function. If this filtering operation is 
applied to the incompressible Navier-Stokes equations one 
obtains 
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Momentum ; 


(3) 


3u. 

1 

3t 


+ 


u . 
3 


9u. 

1 

3x. 

3 



8 

8x . 
3 


T . . + 

13 


a2- 

8x . 8x . 
3 3 


Continuity: 


8u . 
1 

8x. 


0 


(4) 


Here is the subgrid scale (SGS) stress tensor, and 

density, p, is imbedded into the definition of p. The equations 
will be closed once is modeled. In reference 3 for the large- 

eddy simulation of turbulent channel flow, an eddy viscosity 
type model is used for i this model, the eddy viscosity, 
v^, is directly written as a function of , which is the 
strain rate tensor of the large-scale field. Owing to this 
strong coupling between the SGS-stresses and large-scale field, 
this model absorbs too much energy from the large-scale motion 
and inhibits evolution and development of turbulent flow. The 
prescription of an appropriate model which represents the 
correct dynamics of the flow field will be a major part of our 
second-year work. 

The main task performed was the development of a vector- 
ized computer code for efficient numerical integration of 
equations (3) and (4). These equations are written in a form 
which will conserve energy, e.g. 

Momentum; 



8u . 
1 



0 


( 6 ) 
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Here Re is the Reynolds number based on channel centerline 
velocity, U , channel half-width, 6, and kinematic viscosity, v. 

is a function that contains the convective terms, i.e., all 
the SGS terms along and and part of the SGS terms, i.e., all 
For convenience, we have dropped the overbar from the filtered 
variables. We have used a semi-implicit algorithm to numeri- 
cally solve equations (5) and (6) . This procedure employs the 
following finite difference formulae for time -advancement. 

i) For H^: Explicit Adams-Bashf orth 

ii) For 9p/8x.: Implicit Crank -Nicolson 

2 ^ 2 2 2 2 2 

iii) For 3 u^/9x^, 9 u^/9x2, 9 u^/9x2: Implicit Crank— 

Nicolson 

iv) is evaluated at the previous time-level 

The solution procedure assumes the flow to be homogeneous 
along x^^- and X 2 ~directions and makes use of the pseudo-spectral 
technique to calculate space derivatives along these directions 
via fast Fourier transforms (FFT) . Three-point central- 
difference operators are used along the X 2 -direction. Figure 1 
shows a schematic of the flow field geometry. 

Once the governing equations are discretized in time, we 
employ a 2-D FFT along x^ and x^ and transform the equations 
into the wave-number space, k^, k^ . This transformation enables 
one to obtain a set of ordinary differential equations along X 2 
which is solved for every pair of k^ and k^ . 

Hence to solve both imaginary and real parts of the 
dependent variables, two sweeps through the data base are 
required. During these sweeps a block-tridiagonal system must 
be solved for each pair of k^ and k^ along the X 2 -direction. A 
special solution technique is employed for k^ = k^ = 0, for 
which a scalar tridiagonal system is solved separately. 

It should be noted that the solution procedure involves 
three main mathematical operations which had to be vectorized. 
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(1) One-dimensional Fast Fourier Transform: This is used 

to evaluate the derivatives along the x^- and X 2 “directions by 
the pseudo-spectral technique. For this purpose we have imple- 
mented the STAR library subroutine Q4FF0RMS, which is recom- 
mended to be used for a large number of systems, each consisting 
of a large n\amber of independent transforms. 

(2) Two-dimensional Fast Fourier Transform: For the trans- 

formation of the difference form of equations (5) and (6), we 
have employed the STAR library subroutine Q4FFT2DR, which mini- 
mizes array storage by making use of the two-dimensional Fourier 
transform properties of a real field. 


(3) Block-Tridiagonal matrix inversion along X 2 : Note that 

the governing equations can be written in the following vector 
form. 


A - B p-l . C p-i 


= R 


(7) 


In equation (7) F is the solution vector whose components 
are u^^, p; R is the right-hand-side vector that contains the 

convective terms and part of the viscous terms. A, B, C are 
coefficient matrices whose components are functions of the kine- 
matic viscosity, v, and wave numbers k^^ and k^- The solution 
procedure requires the inversion of this block-tridiagonal system 
along X 2 for each k^^ and k^ in two passes, once for real parts of 
u^,U 2 imaginary parts ^ 2 /P and once for imaginary parts of u^,U 2 
real parts U 2 /P. In the vectorization process we have used a 
vectorized subroutine developed by J. Lambiotte (in a private 
communication) which enables the inversion of a large number of 
tridiagonal systems simultaneously decreasing the number of scalar 
operations required to invert each system separately. This 
results in a very substantial reduction in CPU time. In imple- 
menting this subroutine we have found it to be most efficient to 
do the inversion at two x^^ - x^ planes simultaneously only once 
for both passes. This was done by a DO-LOOP such that at each 


7 


index of the DO LOOP a total of 62 tridiagonal systems are 
inverted. 

In addition to these, the rest of the program has been 
written in STAR-FORTRZUsI which makes extensive use of explicit 
vector statements to decrease CPU time. The current program 
requires about 5 seconds of CPU time and 3 CRU per time-step 
on the CYBER-203, including input-output. 


3 . RESULTS 

3.1 Test Cases for the Plane Poiseuille Flow 
(No Suction-Blowing) 

To test the computer code developed during the first year, 
we have performed a number of preliminary computations. The 
results of the cases without suction-blowing are contained in 
this section. 

As the first test we have chosen the calculation of entry 
flow in a channel as a model problem and prescribed the 
following initial conditions at Re = 100 


u, = 


C(1 - X 2 ) + eL^ sin(irx 2 ) cos 


X 1 

4 t \ -ij 

sin 

^3l 

2tT 



L 3 J 


Uo = 


- e(l + cos 'n'X 2 ) sin 


4^ L~ 

sin 

^3l 

iT 



L 3 J 


TT 


U- 


-eL^sin 


47T = — sin(TTx^) cos 
1 ^ 


^3' 

iT 

^3- 


( 8 ) 

(9) 

(10) 


where L^^ = 2tt and L^ = 4 tt/ 3 are the box lengths along the 
homogeneous directions x^ and x^, respectively. The value of 
the constant, e, was taken equal to about 5% of the maximum 
velocity in the channel. 
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Results of this case are presented in figures 2 and 3. 

Figure 2 shows the evolution of the mean velocity profile, <u^>, 
and indicates that with increasing time-steps the parabolic 
velocity profile of the fully developed channel will be 
recovered. The time-history of the turbulent statistics (ref. 1) 
shows that the disturbance energy rapidly decays with time as 
expected. 

During this code verification phase we have performed 
a second calculation by using the two-dimensional Orr-Sommerfeld 
solutions as the initial conditions. The complex frequency, w, 
is an output of the solution procedure and represents the 
eigenmode of the flow for given a, 3, Re and U(x 2 ). For this 
calculation we have selected Re = 1500, a = 1 and 3=0. The 
Orr-Sommerfeld code (ref. 2) gives the distribution of U 2 
and 03, CO = -0.32630 + 0.02821 i, as the solution for this case. 

Since u = li. (y)e^'^^^ , the imaginary part of co determines 

^ ^ —0 0281 
the nature of the solution. For this case, e ‘ , and 

the solution decays in time. 

Using the initial conditions generated as explained above, 
we have performed a series of calcijlations using our code. 

Along the X 2 ~direction we have used a mesh clustering of the 
form 


Yj = ^ tanh[rij tanh ^ (a) ] 

where Yj is the X 2 coordinate in the physical plane and rij is 

the X 2 coordinate in the stretched plane (uniform mesh spacing) . 

In these calculations we have used A == 0.8, which gives 

(Ax)^. 0.035. 

min 

It should be noted that the time-advancement scheme we 
employ is partly explicit and partly implicit. Although in 
light of linear stability analysis, implicit methods are 
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sometimes advocated to be unconditionally stable, time-accuracy 
requires adherence to stability bounds of explicit schemes. In 
all the calculations reported here , we have obeyed both the 
convective stability requirement (the Courant-Friedrichs-Lewy 
conditions) and the diffusive stability requirement. At each 
time-step for the former we used 


C.N. 


AT 



+ 

^2 


^3 

I— 1 

<1 

Ax2 


Ax^ 


max 


< 


0.35 


( 11 ) 


where C.N. is the Courant number. For diffusive stability we 
used 


D 


Re 


AT 


Ax: 


< 0.1 


max 


(12) 


With (AXo) . ~ 0.035 we have found it sufficient to use 

2 'min 

0.02 AT ^ 0.04 to satisfy both of the above stability criteria. 

In choosing the periodicity length (the box length) v/e set 
the smallest wave number allowed in the computational domain 
equal to a (along x^^ and x^) • This corresponds to the longest 
wavelength allowed, i.e.^ the one whose wavelength is equal to 
the periodicity length. 

Results for this case are summarized in figures 4 and 5. 
Figure 4 shows the plot of u^(x 2 ) at two instances in time; the 
solid curve corresponds to the initial disturbance at T = 0 and 
the broken line is at T = 260, which corresponds to a phase 
shift of IT. The shape of the disturbance distribution is well 
preserved and a strong decay is also indicated. Accordingly the 
computed phase angle is equal to 0.291, as opposed to the exact 
solution, 0,3263 (which is the real part of o)) . 
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In figure 5, time-history of plane averaged disturbance- 
velocity components are plotted. These show that the distur- 
bance energy decays almost exponentially (as expected from the 
linear theory) except in the region T < 2.5. The computed 
results decay at a rate about twice as fast as the Orr-Sommerf eld 
solution shown as the solid line. The reasons for this V 7 ill be 
explained later. 

A series of calculations with finite-amplitude disturbances 
were also made for a test case which was run for identical input 
conditions as Run No. 2 of reference 3. These conditions are 
summarized in Table 1. 

Results for this case are contained in figures 6 to 9 . 

Figure 6 shows the time-history of average disturbance inten- 

/N 2 /N 2 

sities, and and indicate a strong damping of these 

quantities. Figure 7 shows a comparison of the present calcula- 
tions with those of reference 3 for the time-history of the 
maximum amplitude of u. The computations of reference 3 
indicate that after an initial transient period, the disturbance 
grows slowly. Our calculations, however, show a strong damping 
of the initial disturbance. This behavior of our computations 
is very similar to the previous case in which the decay rate 
of small amplitude disturbances was found to be much faster 
than indicated by the Orr-Sommerf eld solution. 

/v2 

Figure 8 shows the time-history <U 2 ^ indicates a fast 

growth for this quantity. This lends some support to the 
findings of reference 3, such that an initially two-dimensional 
finite amplitude disturbance may be driven to instability by 
three-dimensional effects caused by round-off errors. In this 
case, the most likely source of round-off errors are the Fourier 
transform operations. 

Finally, in figure 9 a plot of profiles of <u^> is 
presented at three instants in time. The shape of the initial 
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profile is well preserved, and a strong amplitude decay is 
displayed. 

The test cases presented herein have indicated that our 
results are more dissipative than both the Orr-Sommerfeld solu- 
tions and the results of reference 3. It should, however, be 
noted that the method of reference 3 employs Chebyshev poly- 
nomials along the cross-stream direction so that the accuracy 
of spatial derivatives is "infinite order" with respect to 
the second-order finite-difference scheme that is employed in our 
method along the cross-stream direction. Hence, we believe that 
the dissipative nature of our code is caused by the filtering 
effects of the finite-difference scheme. In fact, in compu- 
tations of the complete transition process in which the 
governing equations will be filtered, this dissipative property 
will provide an implicit "filter" along the cross-stream 
direction. The governing equations will be formally filtered 
along the homogeneous directions only. For the case in 
which stability of laminar flow in the linear range is of 
interest, a more accurate method along the cross-stream 
direction seems necessary. This can be obtained either by 
increasing mesh points along this direction or by employing 
polynomial expansions for the flow field variables. 


3.2 Test Case with Suction-Blowing 

In order to implement suction boundary conditions, it must 
be recalled that our solution procedure employs homogeneous 
boundary conditions in the streamwise and the spanwise x^ 
directions. This, of course, requires that the incoming mass 
flow rate must be equal to the outgoing mass flow rate. Hence, 
the continuity equation averaged over x^ - x^ planes gives 


9x2 


0 


(13) 
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from which 


<U 2 > = const. (14) 

Equation (14) states that continuity can be preserved only 
if the average velocity at one wall is equal in both magnitude 
and direction to the average velocity at the other wall. 
Therefore, with this formulation, the only mathematically 
plausible physical problem is to apply suction at one wall and 
blowing at the other wall with equal magnitude. 

The velocity boundary conditions take the form u^^ = 0, 
u^ = 0, and U 2 = f (x^,X 2 ) where f(x^^,X 2 ) can be any periodic 
function. The pressure boundary condition at each wall is 
obtained from the U 2 momentum equation evaluated at the first 
interior point away from the wall. An exact solution for 
plane-Poiseuille flow with suction-blowing boundary conditions 
can be obtained by integrating the equations of motion for 
steady, incompressible flow in a two-dimensional channel. This 
gives , 


/ 

Re V 1+x^) 

/N 1 

2Re V 

\ 


o c 2 

2e - e 


/N 

— ^ 'V 


(15) 


Re V. 


2ReV„ 

^^2 

e - 1 

^ > 

> 


Here the caret refers to a nondimensional quantity, and Re 
is the Reynolds number based on channel half-width and channel 
centerline velocity. 

It should be noted that equation (15) becomes interdeter- 

/s 

minate for V„ = 0 which simulates flow in a channel without 
c 

suction-blowing . However , 
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(16) 


2 

lim u, = (l-x„) 
V ^ 0 1 2 

c 


Hence as V^-^0 expression (16) recovers the parabolic velocity 
profile of the plane Poiseuille flow. 

The suction-blowing boundary conditions were tested by 
performing calculations at Re = 100, = 0.05, where is the 

suction velocity nondimensionalized by channel centerline 

velocity. The following initial conditions were used. 

u^ = C(1 - x^) 



Note that these conditions simulate a fully developed 
channel flow which is subjected to suction at one wall and 
blowing at the other wall with equal velocity, V^. Calculations 
were performed with T = 0.006 (as in section 2.2) and results 
are presented in figure 10 at T = 0.3 which corresponds to 
50 time-steps. In this figure the exact solution is also 
presented, and a very good agreement between the two solu- 
tions is displayed. As expected, in both cases the average 
velocity distribution is skewed and the location of maximiim 
velocity has moved towards the suction-wall. The velocity 
gradients (in both cases) are considerably large at the suction- 
wall than at the blowing wall. 
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4. 


STUDY OF SUBGRID SCALE MODEL AND 
LENGTH SCALE 

During this contract period we have made a preliminary 
study of a subgrid scale (SGS) turbulence model to be imple- 
mented into the code for the full simulation of transition. It 
should be noted that the type of SGS models in which the eddy 
viscosity is written as a function of the strain-rate of the 
large-scale motion extracts too much energy and inhibits the 
flow development. Considering the strong coupling between the 
filtered and unfiltered fields due to this type of SGS model, 
the use of a SGS model which relaxes this direct coupling will 
likely improve the calculation of transition. Moreover, it is 
desirable that a SGS model should be able to "turn on" itself 
due to some evolving property of the flow field which is 
characteristic of turbulence. Such an evolving quantity can 
be represented conveniently in terms of the kinetic energy of 
small-scale motion. 

In this section, we formulate a SGS model which is derived 
from the transport equation for the kinetic energy, E, of the 
small-scale motion. We consider the equation 


3E 

9t 


+ 



9E 


convection 


i!|- 

f9u.' 

1 

2 

f9u/ 

1 

2 ~ 

1 

'9u^ 9u.' 



+ 



2 ^ij 

9x7 9x. 

3 iJ 


■ ■ S/- ■>-■■■■ ■ 

dissipation (s) production (P) 


2-2 
u, u . u, u . 
k 1 k 1 


- 2 - 




9x, 


_ -U.U, U. + U.U, + 

2 iki Ik PP 


(17) 


velocity and pressure diffusion 
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The turbulent channel flow is an equilibrium flow in which 
the contribution of the transport terms to the energy balance 
is much less than the other terms. Also since the effect of the 
SGS terms is to provide a physically plausible "sink" term for 
energy lost from the large-scale motion, a detailed and time- 
accurate solution of the kinetic energy equation is likely not 
necessary. In this work, therefore, we intend to use the 
simplest form of equation (17) with temporal and transport 
terms neglected. This gives: 

Dissipation = Production 

With this simplification the only higher-order term to be modeled 
is the dissipation. To this end, we introduce the following 
relation 


e 



[1 + C^/Re] 


(18) 


where 


Re 




and Cj^ and are constants to be determined. 

Note that formula (18) has the correct forra near the wall, 
i.e. , where Re << 1, and outside the wall region, i.e., where 
Re >> 1. 

Another assumption is required to relate (which is the 

quantity we need to calculate) to E (which is the quantity 
equation (17) calculates) . For the constitutive relation 
between E and t . . we propose to employ a Prandtl-Kolmogorov 


16 



type foonnula such that 


T . . = - 2 v_S . . 
ID T 13 


and the eddy viscosity is given as 


Vt = C2LE 


1/2 


(19) 


( 20 ) 


where C 2 is a constant to be determined from physical consider- 
ations. Finally a modification for the presence of the wall is 
needed to suppress the eddy viscosity very near the wall. This 
can be obtained by multiplying by a function f ^ , which can 
be defined as 


fy = [1 - exp (-0.014 Re)] 


( 21 ) 


It should be noted that the use of the kinetic energy 
equation requires the definition of a length scale. To this 
end we intend to implement the following formula for L, 

L = ^min(A^,^,') min(A 2 ,il') min(h 2 ,il')j (22) 


where i' is Prandtl's mixing length (ref. l). 


0.1 

Y > 

0.1/k 

KV 

y < 

0.1/k 


(23) 


In equation (22) A^ and A^ are filter widths in the x^^- and x^- 
directions, h 2 is the grid size in the X 2 ~direction, y^ is the 
distance to the nearest wall and k is the von Karman constant 
equal to 0 . 4 . 
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A final assessment of the model and final evaluation of 
model constants will be made . 


5. CONCLUSIONS 

A vectorized, semi-implicit code for simulating transition 
and laminar flow control by suction in plane Poiseuille flow 
is developed. The code takes about 5 sec. of CPU time and 3 CRU 
per time-step on the CYBER 203 for a grid of 32x33x32 mesh 
points . 

The code has been tested with various boundary conditions 
(no-slip, suction-blowing) and initial conditions [(a) arbitrary 
periodic, (b) Orr-Sommerf eld eigensolutions] and has been found 
to be stable and fast-converging. Comparisons of the present 
results with the linear theory have shown that the present 
results are more dissipative than both the Orr-Sommerf eld 
solutions (the linear theory) and the results of reference 3. 

We attribute this to the "filtering effect" of finite-difference 
operators which are employed along the nonhomogeneous direction. 
For the full simulation process this "dissipative" property will 
provide an implicit filter along the cross-stream direction. 

The governing equations will be formally filtered along the 
homogeneous directions only. For the case in which stability 
of laminar flow in the linear rage is of interest, a more 
accurate method along the cross-stream direction seems necessary. 
This can be obtained either by increasing mesh points along this 
direction or by employing polynomial expansions for the flow 
field variables. 
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Table 1.- Summary of Conditions for Test Case 
with Finite-Amplitude Disturbances 


Reynolds number 

2935 

x^-wave number {a = 2II/x) 

1.3231 

X 2 -wave number 

0 

Reo) 

0.4353 

Imco 

-0.01819 

Initial amplitude 
(x^-velocity) 

.1175 

Spatial resolution 

32x33x32 

Time step, AT 

0.02 

Final time, T 

40 
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Tigure 1.- Flow field geometry, the computational box. 
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Figure 2,- Entry flow in a plane channel with 3-D disturbances. 




Figure 3.- Time-history of disturbance energy 



Figure '4.- A plot of the profile of U 2 ^(x 2 ) , for the small 

disturbance case; , initial disturbance, — — , at 

t = 10.8, which corresponds to a computed phase 

angle of tt. 
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Figure 5 
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- Tune -history of average disturbance intensities O — * 
A — , < u* >; , Or r- Sommer f eld solution. 






Figure 6.- Time-history of average disturbance 
intensities; finite amplitude case. 











Figure 10.- Laminar flow in channel with suction, 
V /u = 0.01, Re = 100. 
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